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ABSTRACT 

The massive red supergiant W26 in Westerlund 1 is one of a growing number of red supergiants shown to have winds that are 
ionized from the outside in. The fate of this dense wind material is important for models of second generation star formation in 
massive star clusters. Mackey et al. (2014) showed that external photoionization can stall the wind of red supergiants and accumulate 
mass in a dense static shell. We use spherically symmetric radiation-hydrodynamic simulations of an externally photoionized wind 
to predict the brightness distribution of Ha and [Nil] emission arising from photoionized winds both with and without a dense 
shell. We analyse spectra of the Ha and [NII] emission lines in the circumstellar environment around W26 and compare them with 
simulations to investigate whether W26 has a wind that is confined by external photoionization. Simulations of slow winds that 
are decelerated into a dense shell show strongly limb-brightened line emission, with line radial velocities that are independent of 
the wind speed. Faster winds (> 22 kms^*) do not form a dense shell, have less limb-brightening, and the line radial velocity is 
a good tracer of the wind speed. The brightness of the [NII] and Ha lines as a function of distance from W26 agrees reasonably 
well with observations when only the line flux is considered. The radial velocity of the simulated winds disagrees with observations, 
however: the brightest observed emission is blueshifted by 25 km s“* relative to the radial velocity of the star, whereas a spherically 
symmetric wind has the brightest emission at zero radial velocity because of limb brightening. Our results show that the bright nebula 
surrounding W26 must be asymmetric, and we suggest that it is confined by external ram pressure from the extreme wind of the nearby 
supergiant W9. We obtain a lower limit on the nitrogen abundance within the nebula of 2.35 times solar. The line ratio strongly favours 
photoionization over shock ionization, and so even if the observed nebula is pressure confined there should still be an ionization front 
and a photoionization-confined shell closer to the star that is not resolved by the current observations, which could be tested with 
better spectral resolution and spatial coverage. 

Key words. Stars: winds, outflows - Hydrodynamics - radiative transfer - methods: numerical - individual objects: Westerlund 1 - 
individual objects: W26 


1. Introduction 


Most stars form in clusters ( |Lada & Lada|2003 1 and, in particu¬ 
lar, massive stars are predominantly found in clusters and asso¬ 
ciations ( Gies|1987l l. Furthermore, many O stars not currently in 
a cluster have been identified as runaway stars (e.g. Gvaramadze 
et al. |2012| l. Feedback from massive stars in clusters is impor- 
tant for the enrichment and evolution of the interstellar medium 
(ISM) of gas-rich galaxies ( |Krumholz et al.|2014) l, and it is a cru¬ 
cial ingredient in models of multiple stellar populations in glob¬ 
ular clusters (e.g. |Gratton et al.|2012| l. Globular clusters formed 
early in cosmic time, during the epoch of first galaxy formation, 
and their multiple populations (e.g. Piotto et ar]|2007 1 demon¬ 
strate that somehow gas was able to accumulate, cool, and be¬ 


come gravitationally unstable within the star cluster (D’Ercole 
et al.|2008) Palous et al. 20141. It is not at all clear how this hap 


pens; a number of theories and their shortcomings are discussed 
in depth in Bastian et al. ( 2015[ ). It is fair to say that there is 
no consensus on how globular clusters formed. Because of this, 
researchers have begun to look at present-day massive star clus¬ 
ters to try to gain some insight (e.g. Bastian et al.|2013l Bastian 
|& Strader|2014|l, with the result that little evidence of cold gas 


and dust has been found in the intracluster medium of nearby 
massive star clusters. 

Models for multiple populations that involve mass loss from 
massive stars (as opposed to asymptotic giant branch stars, e.g. 


P’Ercole et al. 

2008)1 all concentrate on winds from hot stars 

(e.g. 

Canto et a 

2000 IWiinsch et al.|2008l|de Mink et al.|2009l 

[Krause et al.||20 

13 1 . The difficulty these models face is how to 
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cool down the hot shocked wind to temperatures low enough 
to allow star formation ( |Palous et al.|2014) l. Surprisingly, these 
models do not explicitly consider red supergiants, even though 
Galactic stars with initial mass M < 30 Mq lose more mass 
as red supergiants than they do as main sequence stars ( |Chiosi| 
|& Maeder||1986| l. Their winds are already cold (temperature 
ss 100 - 500 K; e.g. Schuster et al.|200^ Le Bertre et al.|20T^ , 
neutral, dusty, and dense, and so if this wind material can remain 
in a star cluster then it is the obvious candidate to host second 
generation star formation. Furthermore, evolutionary models of 
very metal-poor stars with initial masses above 100 Mq show 
that they may evolve to red supergiants already on the main se¬ 
quence (Marigo et al. 2003 Szecsi et al. 2015[ Kohler et al. 
|2015[r%on et al.||2012 l. For this reason we consider it impor¬ 
tant to study mass loss from red supergiants in star clusters, to 
constrain whether their wind material can remain cold and dense 
in such a harsh environment. 
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A growing number of red supergiants are found to be sur¬ 
rounded by nebulae with emission lines from ionized gas, show¬ 
ing that their winds are photoionized by external radiation fields. 
This was first proposed for NML Cyg ([Morris & Jura 1983|l, then 


and more recently W26 ( 

Wright et al.[[2014[l and IRC -10414 | 

([Gvaramadze et al.[ 

2014 

[Meyer et al. 2014a[l. WOH G64 in if 4 x 10 


the Large Magellanic Cloud also shows similar nebular emis¬ 
sion lines (Levesque et al.|2009 1. Mackey et al. ( |2014[ l proposed 


that Betelgeuse’s wind is also externally photoionized, to explain 


the presence of a static neutral shell around the star (Le Bertre 
et al. [20 12 ) 1 . We expect on this basis that many red supergiants 


Westerlund 1 (|Westerlund| 1961 1 


_ - 

massive young star cluster in our Galaxy ( |Brandner et al.|20d8| l, 
and its age (4+1 Myr) is such that it contains a population of 
extreme supergiant stars ( jClark et al.||2005) . Only at this age 
do star clusters contain the hottest (Wolf-Rayet stars) and cold¬ 
est (red supergiants) stages of massive star evolution simultane¬ 
ously. The red supergiant W26 (RA: 16 47 05.403; DEC; -45 
50 36.76; SIMBAD identifier: Cl * Westerlund 1 W 26) is sur¬ 
rounded by a bright radio nebula ([Clark et al.|1998[ 




FORS spectra - 

MARCS model- 


6000 


have photoionized winds. [Mackey et ar] ( |2014| ) showed that the 
winds of red supergiants could be decelerated and accumulated 
into a dense shell, if the wind is exposed to an external ionizing 
radiation field from all sides. This could have important conse¬ 
quences for the fate of the mass lost during the red supergiant 
phase, and also for any subsequent evolutionary phases and the 
observational characteristics of the star’s eventual supernova ex¬ 
plosion ( [Mackey et al. 2014 1 . If this mechanism is effective, then 
the mass lost by red supergiants could remain within the parent 
star cluster for much longer than if the wind is freely expand¬ 
ing, increasing the likelihood that secondary star formation could 
happen. 


7000 8000 

wavelength [A] 


9000 


198711 is probably the most 


Doughert^ 


|1998[ l and 


[et al.[[20T0] l, also seen in mid in frared ([Clark et'af 
in narrow-band optical emission ( Wright et al.|2014 i. The wind 
of W26 is so bright in typical tracers of photoionized gas (Ha 
and [Nil] forbidden-line emission) that a number of authors 
have suggested that it could be externally photoionized by ra¬ 
diation from the O-type and Wolf-Rayet (WR) stars in the clus- 
ter ([Dougherty et al.[[2010| [Wright et al.[[2oi4[ Mackey et al.[ 
2014[ l. The main alternative explanation for the ionized nebula 
around W26 is that it is hydrodynamically confined by the pres¬ 
sure of the intracluster gas. This could be ram pressure from stel¬ 
lar winds or a global cluster wind ( [Dougherty et al.[|2010 i, or 
thermal pressure from thermalised shocked wind bubble and re¬ 
cent supernovae that must have occurred in Westerlund 1 ( [Muno[ 
et al.|200^ Kavanagh et al.|20TT[ l. 


Fig. 1. Comparison between observed FORS spectra for W26 (black 
solid line) and MARCS model fluxes (dashed red line) for a = 3600 
K supergiant. The wavelength gap where no spectrum was obtained 
arises from a gap between the two grisms. 


To make progress, we take the existing spectroscopic data on 
W26 and its circumstellar environment, and compare this to the 
predictions of numerical simulations of an externally photoion¬ 
ized wind, following the method of Mackey et al. ( 2014[ l. The 
main observables are the radial velocity of the emitting gas, the 
intensity of this emission, and the line ratio, as a function of po¬ 
sition along the slit in the Ha and [NII] lines (in this paper we 
focus entirely on the brighter [Nil] line at 6584 A). As shown 
below, these data already provide contraints that are difficult to 
satisfy with a spherically symmetric photoionized wind. 

The observational data, reduction methods and results are de¬ 
scribed in Sect. [2] The results of numerical simulations are de¬ 
scribed and compared with observations in Sect. We discuss 
our results in Sect.[4|and conclude in Sect. [5] 


External photoionization and hydrodynamic confinement 
must both be present, because the combined winds of the evolved 
stars will create a high pressure environment ( Muno et al.|2006] l 
and there are many hot stars in the cluster that emit ionizing pho¬ 
tons, e.g., 24 WR stars ( [Crowther et al.[[2006[ ). The challenge 
then is rather to understand which process is producing the ion¬ 
ized nebula in W26, and by consequence in the winds of red 
supergiants in other star clusters. The radius at which a stellar 
wind becomes ionized is well defined for a specified external ra¬ 
diation field ( [Morris & Jura|1983] l, and the hydrodynamic con¬ 
finement radius is set by the balance between the ram pressure 
of the stellar wind and the thermal/ram pressure of the environ¬ 
ment. Unfortunately there are no good estimates of the mass-loss 
rate or wind velocity of W26, and we do not know the physical 
distance separating it from the other stars (such as the supergiant 
sgB[e] star, W9) that are nearby in projection. 


2. Dataset and analysis processing 

W26 was observed with the EOcal Reducer and low dispersion 
Spectrograph (EORS 2) (Appenzeller & Rupprecht|[1992[ Ap-[ 
[penzeller et al. [ 1 998| ) attached to the Cassegrain focus of the ESO 
Very Large Telescope of the Paranal Observatory, Chile. The ob¬ 
servations were carried out on the 16th of April changed 2011 
under the ESO program ID 087.D-0673(A). The data were ob¬ 
tained in long-slit mode using the 1200R and 1028Z grisms and, 
a slit size of 0.3 arcseconds to obtain the highest resolution avail¬ 
able. This instrumental configuration offers a wavelength range 
between ~ 5700 - 9300 A and a resolving power of R ~ 7000 
( [Negueruela et al.|2010[ [Clark et al.|2013[l. Th e slit location and 
orientation is shown in fig. (2) of [Wright et al. ( 2014[ l, where the 
data were first presented. 

The data reduction was performed employing IRAl|^ ( |Tody| 
1993[l and IDL routines following the standard steps and tasks 


for long-slit spectral data reduction (consisting of bias subtrac¬ 
tion, fiat-fielding, spectral extraction with and without back¬ 
ground subtraction, and wavelength calibration) and analysis. 


* Image Reduction and Analysis Facility (IRAF - 
http://iraf.noao.edu/1 is distributed by the National Optical 
Astronomy Observatory, which is operated by the Association of 
Universities for Research in Astronomy (AURA) under cooperative 
agreement with the National Science Foundation. 


Article number, page 2 of 19 




















































































































Mackey et al.: Cold gas in hot star clusters: the wind of W26 




-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 

Distance to W26 (arcmin) 


Fig. 2. Spectral observations of W26 and neighbouring stars, (a) radial 
velocity of the Ha (black dots) and [NII] (blue triangles) lines along the 
slit in the W26 co-moving reference system (-32kms“'), with posi¬ 

tion measured relative to W26. The directions along the slit are also in¬ 
dicated: WNW for Westnorthwest and ESE for Eastsoutheast. (b) Line 
brightness (intensity) in the two spectral lines along the slit (arbitrary 
units), (c) Fitted Gaussian linewidth, <t in the two spectral lines along 
the slit. The three stars included in the slit, and no further away than 0.4 
arcmin from W26, are labelled according to ([Negueruela et al.|2010^. 


2.1. W26 Stellar parameters 

We used the background-subtracted, wavelength-calibrated 
FORS2 spectra of W26 to roughly estimate the effective tem¬ 
perature (Teff)- We did this by comparing the observed spectrum 
with synt hetic fluxes calculated with marcs stellar atmosphere 
model^ ( [Gustafsson et al.||2008| with spherical symmetry and 
solar metallicity ( [Grevesse & Sauval||1998| spanning between 
3300K and 4500K in Tefr, in steps of 100K, and -0.5 and -i-l.O 
in log g, in steps of 0.5 dex. 

Given the very low resolution of the FORS2 spectrum, hence 
the heavy line blending, there was no continuum available that 
could be used to properly normalise the observed spectrum. For 
this reason, we normalised the observed spectrum using each 
model flux at a time as a reference before determining the total 
deviation between observed and model spectrum. We obtained a 
best htting effective temperature of T^ff - 3600 K, which corre- 

^ http://marcs.astro.uu.se/ 


sponds to a spectral type of about M2.5 (|van Belle et al.|1999 , 


in agreement with previous spectral type determinations (Clark 
et al.|2010)l; a comparison of the model with the FORS spectrum 


IS shown in Fig. [T] The formal uncertainty in Tefr is 200 K, but 
systematic uncertainties in the modelling are likely to be larger. 
It is important to recall here that W26 is a known variable star 
with a spectral type varying between M2 and M5 (IClark et al.l 
[ 2 ^ . 


2.2. Radial velocity and Intensity distribution maps 

The kinematic study presented here is based on the the radial ve¬ 
locity distributions of Ha 6562.79A and [Nil] 6583.45A along 
the spatial direction of the FORS 2 long-slit spectra. This analy¬ 
sis was done with the 1200R grism data since these two transi¬ 
tions can only be reached in its wavelength coverage. To explore 
the entire spatial information available in the spectra, each row 
of the science images was sliced and individually calibrated by 
the equivalent row in the calibration arc-lamp spectrum. The tar¬ 
geted lines were htted by a Gaussian function per row, identify¬ 
ing the position of the emission peaks and their amplitudes with 
respect to the spectral continuum. The velocities are presented in 
the heliocentric reference system. Previous studies reported the 
possible presence of small instrumental flexures that migh t lead 
to shifts in the wavelength calibration as large as 5kms * (Bag- 
|nulo et al.|2013| l. We corrected for such a possible instrumental 
bias by further calibrating the spectrum using two telluric sky 
emission lines at 5889.951 and 6300.304 A. 

Fig. 0 displays the radial velocity (vios, panel a) and line 
brightness (b) spatial distributions of Ha and [Nil] along the 
FORS 2 slit up to a maximum distance of 0.4 arcminutes from 
W26. The radial velocity was shifted to the W26 radial veloc¬ 
ity co-moving frame, -32kms“* (Negueruela, private com¬ 
munication), obtained using data from the AAOmega spectro¬ 
graph on the Anglo-Australian Telescope, and reduced using the 
techniques described in |Marco et al. ( |2014| l, with a precision of 
2 - 4kms“'. We derived uncertainties of a; 4 and 7kms“' for 
Ha and [Nil] radial velocities, respectively. 

We also measured the linewidths of the Ha and [NII] lines 
along the slit, shown in panel (c) of Fig. and from the Gaus¬ 
sian hts we obtained cr a; 30 km s ' for both lines. For Ha most 
of the measurements are in the range 30-40 km s ', whereas for 
[Nil] they are mostly in the range 20-40 kms“*. This is consis¬ 
tent with the instrumental resolution of FORS 2. The only ex¬ 
ception is near W9, where the very bright Ha line has a width 
cr ^ 55 kms *. There is an apparent trend that cr decreases 
from left to right across the nebula around W26, or with increas¬ 
ing distance from W9. It is unclear if this is signihcant because 
we are so close to the instrumental resolution; higher resolution 
spectra could conhrm this. Such a trend might be expected if the 
wind from W9 were directly colliding with the wind of W26 (see 
section]^. 

For photoionized gas at T a; 10"' K, Ha has a thermal 
linewidth cr =» 10 kms *, and the [Nil] line is signihcantly nar¬ 
rower. We therefore cannot hope to resolve any substructure in 
the nebular emission, for example blue- and redshifted compo¬ 
nents of the wind of W26. This means that the only informa¬ 
tion we can reliably obtain from the spectral lines are the radial 
velocity of the peak and the line brightness. It should be noted, 
however, that the observed line could be a blend of multiple unre¬ 
solved lines with differing radial velocities. Furthermore, the line 
brightness is uncalibrated which makes absolute flux/luminosity 
measurements impossible, and so we are limited to comparing 
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the relative line brightness at different locations along the slit. 
Spectroscopy with significantly higher spectral resolution, able 
to resolve multiple line components in photoionized gas, would 
provide much stronger constraints on models. 

W9 has an enormous Ha luminosity (|Westerlund| 1987 1 with 


a measured equivalent width W - -520 + 30 A (Clark et al. 


|2005|l, which they point out is unprecedente d for an early type . 


Clark et al. 


2013| l. 


The three data points closest to W26 have significantly larger 
uncertainty than the surrounding points, because the continuum 
flux from W26 itself fills in the line at radial velocities near zero. 
This biases the radial velocity measurement of the line to the 
blue, and so these three points blueshifted by 30-50 km s ' are 
probably spurious (the line intensities in these points are simi¬ 
larly more uncertain than neigbouring points). Aside from these 
points, the radial velocity of the nebula around W26 is almost 
uniform at a; -20 km s“'. 

In the low-density limit, assuming that all N is in the form 
of and that H is fully ionized, the emissivity ratio of [NII] to 
Ha can be expressed as ( |Dopita|1973| l 

iH. 0.116 VfexpM^ ,1) 

jHa \ T I 6.8 X 10-5 

where we assume T » 350 K, and where /(N) is the number 
fraction of N (relative to H) for which 6.8 x 10 ^ is the solar 
photosphere abundance ( |Asplund et al.|2009) . In our models we 
assume /(N) = 2x10 appropriate for nitrogen-enriched winds 
of rotating red supergiants ( |Brott et ar]|201 1 1, leading to a pre¬ 
factor of 0.341 instead of 0.116. For typical H ii region tempera¬ 
tures (7000-10000 K) the ratio is between 0.5 and 1.5 for solar 
abundance, but always > 1 for the enhanced N abundance. 

The low density limit applies (neglecting collisional de¬ 
excitation) because in our simulations the electron number den¬ 
sity, We, is always less than the critical de nsity for the 6584A line 

table 
5x105 


of [NII], which is W p^rit = 8.6x10"' cm ^ ( Osterbrock 1989 


3.11). For reference Dougherty et al. (2010 1 estimate n, 


cm-5 from radio observations, if one assumes a constant density 
in the emitting region. Comparing this estimate to the critical 
density, we expect a reduction in [Nil] emission of < 10% be¬ 
cause of collisional de-excitation. Taking into account the pos¬ 
sibility of doubly ionized N and of collisional de-excitation, any 
constraints we obtain on the nitrogen abundance are lower limits. 


3. Simulations of an externally photoionized wind 


emission line star (or W - -640 + 40 A; 

The Ha line is much brighter than the [Nil] emission lines. 
The diffuse emission in the ISM far away from W9 and W26 
also has stronger Ha than [NII] emission. In contrast, the neb¬ 
ula around W26 has [Nil] emission significantly brighter than 
Ha, as noted by Wright et al. ( 2014| l. This is not generally 
obtained in shock models with solar abundances, which pre¬ 
dict intensity ratios [NII]/Ha < 1 (e.g. |Raymond|[T979 1, and 
so can be a signature of enhanced nitrogen abundance. Stellar 
evolution models predict that the nitrogen abundance in red su¬ 
pergiant winds should be enhanced by a factor of 3-5 ( |Brott| 
|et al.|[20TT] l, and this has been measured in circumstellar neb¬ 
ulae around a number of evolved massive stars including the 
runaway red supergiant IRC-10414 ( [Gvaramadze et aL||2014[ 
Meyer et al.||2014a| and references therein). Alternatively, ion¬ 


ization fronts can show enhanced [NII] emission under certain 
conditions ([NII]/Ha ~ 1 - 3), even for solar abundances ( |Hen-| 
ney et al.|2005| l. 


effects, in particular decelerating the wind into a dense and mas¬ 
sive shell that could contain up to 35% of the mass lost. This 
completely changes the density and velocity structure of the 
wind, and hence its emission properties. The wind structure was 
calculated by Mackey et al. ( |2014| l using analytic and numeri¬ 
cal methods, in both cases assuming an isothermal equation of 
state (with temperatures of T = 10^ and 10'' K in the neutral 
and ionized parts of the wind, respectively). This allowed them 
to predict 21 cm emission from the neutral shell around Betel- 
geuse, but this is a special case because it is so close to Earth that 
the wind can be spatially resolved at 21 cm. 

Winds around other red supergiants have been detected (and 
spatially resolved) in optical nebular lines ( [Wright et al.||2014t 
[Gyaramadze et al.|2014| l. We wish to model these lines but they 
are very senitive to temperature (see Eq. [^, so it is important to 
calculate the gas temperature self-consistently with heating and 
cooling functions, particularly near the ionization front where 
most of the emission is concentrated (e.g. Henney et al.|2005] l. 
The method we use to do this is described below. 

Our models provide the first predictions for nebular emis¬ 
sion from ionization fronts in red supergiant winds, including 
both radiation-hydrodynamics and postprocessing with radiative 
transfer. The wind properties of red supergiants are often poorly 
constrained from both theory and observations (see e.g. the re¬ 
view by jMauron & Josselin|201l} . Optical spectroscopy of neb¬ 
ular lines offers a new avenue to measure the wind velocity and 
density (and hence mass-loss rate, M), through line intensities 
and radial velocities. This is of crucial importance for constrain¬ 
ing the evolution of evolved stars. 


3.1. Input physics for simulations 


We use the pion code ( |Mackey|2012| l in spherical symmetry with 
a wind expanding from the origin (r = 0) and an ionizing radia¬ 
tion (extreme ultraviolet, EUV, with hv >13.6 eV) source at ra¬ 
dius r = oo. The method largely follows previous calculations of 
photoionization-confined shells around red supergiants ( jMackeyj 
et al. |2014|l, except that we use non-equilibrium heating and 


cooling of the gas ( [Henney et al.|2009l l, non-equilibrium ioniza 
tion of hydrogen, and multi-frequency EUV radiation. This cal¬ 
culates gas temperature more realistically than in [Mackey et al. 


P014[l. The im plementation is almost exactly as described in 
[Mackey et al.| ( [2015[ ), where the scheme was used to describe 
expanding Hii regions around O stars; we refer the reader to 
section 2 of this work for details. The rate equation of hydrogen 
ionization is solved including photoionization, radiative recom¬ 
bination, and collisional (electron impact) ionization. We make 
the approximation that helium is singly ionized whenever hydro¬ 
gen is, so the electron fraction is 1.1 times the H^ fraction (for 
helium number fraction nHe/«H = 0.1). We consider a blackbody 
radiation spectrum with radiation temperature = 4 x lO' K, 
in the energy range 13.6-54.4 eV, and also far-ultraviolet (EUV) 
heating radiation following Henney et al. ( 2009[ l with equal pho¬ 
ton flux in EUV and EUV photons. This algorithm for thermal 
and ionization evolution captures the strong cooling in the neu¬ 
tral gas to tens or hundreds of Kelvin and the photoheating to 
T ai 7000 - 10000 K in photoionized gas, including spectral 


Mackey et al. (20141 showed that external photoionization of 


hardening (and extra heating) near the ionization front (e.g. Hen- 
ney et al.|2005[ l, which is important because [NII] and Ha have 


the winds of red supergiants can have strong hydrodynamical different temperature scaling (Eq.[^. 
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3.2. Wind properties for W26 

There is considerable uncertainty in the wind properties of W26 
(e.g. |Wright et al.|2014| l, so we have run simulations for two dif¬ 
ferent mass-loss rates, M -2x 10“^ and 10^“* Mq yr“*, and four 
different wind speeds, Vm = 15, 20, 25, and 30 km s '. Accord¬ 
ing to the |Mauron & Josselin|(20TT]) ref ormulation of the mass- 
loss rate of|Nieuwenhuijzen & de Jag^ (1990 1 , we can estimate 


M = 5.6 X 10 


\105Lo/ 13500 K/ 


Using Tpff = 3600 K (estimated above) and L = 3.5 x 10^ Lq 
( W right et al.||2014 i we obtain M - 4.2 x 10“^ Mgyr''. This 
is only a very approximate estimate, however, as individual red 
supergiants often have mass-loss rates a factor of 3 larger or 
smaller than empirical relations ( Mauron & Josselin|201 l| l. Fur¬ 
thermore, [Mauron & Josselin] ( |2011| l show that measured mass- 
loss rates for the same object often differ by a factor of a few 
depending on the technique used. 

|Fok et aL] ( 2012| l compare near- and mid-infrared photome¬ 
try for W26 with models and estimate M - 2.2 x 10“^ Mq yr“', 
although their model requires W26 to have a total luminos¬ 
ity (10^"^ Lq) that is much larger than the above-quoted esti¬ 
mates. [Dou^ertyjjet^^ (2010 1 study the nebula around W26 
with radio observations, but did not obtain a measurement of 
the mass-loss rate. They estimate that the nebular mass of W26 
is comparable to that surrounding VY CMa ( [Smith et al.|2001|l. 
VY CMa has a time-averaged M ^ A x 10“^ Moyr^^ (Decin 


jet al.||2006| [Mauron & Josselm||201 1|| from gas emission mea 


suremenls, but dust emission gives much larg 

er estimates: M - 

(1.6 - 4) X 10-4 Moyr-' ( 

Danchi et al. 

1994| 

Smith et al. 12009 


cause themass loss is variable. On the other hand, Dougherf}^ 
et al. (20101 estimate M -2x 10“^ Mq yr“' for W237, another 


(presumably coeval) red supergiant in Westerlund 1, and then 
note that it has a similar nebular mass to W26. These results 
show that the two mass-loss rates that we use are reasonable 
estimates for W26, although they are not hard lower or upper 
limits. 

The distance from W26 of the edge of the Ha emission is 
about 0.02-0.04 arcmin on the Westnorthwest side (hereafter 
WNW) and 0.04-0.06 arcmin on the Eastsoutheast side (here¬ 
after ESE); see Eig.[^ We adjust the ionizing flux to set up simu¬ 
lations that have a photoionized nebula at a separation from W26 
that is appropriate for each side. Westerlund 1 contains a num¬ 
ber of Wolf-Rayet and O stars with ionizing photon luminosities 


<2o 


1049 
( 10 ' 


s *, providing ionizing fluxes at W 26 on the order of 
' - lO'^) cm“^ s“* (see section 4.6 1 , so the fluxes that 


we use are about what we expect for the environment of W26. 


3.3. Simuiation parameters 

We have 16 simulations for two mass-loss rates, two positions 
for the ionization front, and four wind velocities. The identifier 
(ID) of each simulation and the parameters used are listed in 
Table In the main text we present only the simulations with 
M - 2 X 10“^ Moyr ' and ionizing fluxes chosen for the ESE 
side of the nebula (i.e. only Voo is varied). We focus on the ESE 
side of the nebula because it extends further from the star and 
hence the data are less contaminated by the central star. The re¬ 
sults from the simulations with M = 10“^ Mq yr“* are very sim¬ 
ilar to those with M - 2 x 10“^ Moyr^' except that the wind 
in the former case is five times denser than the latter, and so the 


Table 1. Parameters used for the 16 spherically symmetric radiation- 
hydrodynamics simulations. In the identifier (ID), M5/M4 refer to the 
mass-loss rates, V15/V20/V25W30 refers to the wind velocity, and W/E 
denote WNW side or ESE side of W26 (i.e. when the ionization front 
is closer to or further from W26). The ionizing fluxes have been set to 
give ionization fronts that match the observations as closely as possible 
in terms of Ha brightness as a function of distance from the star. 


ID 

Voo (kms ') 

M(Moyr'') 

Fy (cm ^ s ') 

M5V15W 

15 

2 X 10-" 

7.87 X 10'" 

M5V20W 

20 

2 X 10-5 

1.29 X 10" 

M5V25W 

25 

2 X 10-5 

8.29 X lO'" 

M5V30W 

30 

2 X 10-5 

5.83 X 10'° 

M4V15W 

15 

1 X 10-4 

1.70 X 10'’ 

M4V20W 

20 

1 X 10-4 

2.34 X 10'’ 

M4V25W 

25 

1 X 10-4 

1.44 X 10'’ 

M4V30W 

30 

1 X 10-4 

9.95 X 10" 

M5V15E 

15 

2 X 10-5 

1.70 X 10'° 

M5V20E 

20 

2 X 10-5 

2.78 X 10'° 

M5V25E 

25 

2 X 10-5 

1.79 X 10'° 

M5V30E 

30 

2 X 10-5 

1.26 X 10'° 

M4V15E 

15 

1 X 10-4 

3.67 X 10" 

M4V20E 

20 

1 X 10-4 

5.05 X 10" 

M4V25E 

25 

1 X 10-4 

3.10 X 10" 

M4V30E 

30 

1 X 10-4 

2.15 X 10" 


line emission is approximately 25 times brighter. It is not possi¬ 
ble to calibrate the data in absolute flux, and so we do not know 
the absolute line brightness and cannot constrain the wind den¬ 
sity (but see Sect.|^. Eor completeness, the other simulations are 
presented in Appendices [A| and [B] 

All simulations use 10 240 uniformly spaced grid zones from 
r = 9.5 X 10 '^ cm 0.0031 pc) to 6.172 x 10 '’ cm ( 0.2 pc). 
This is more than sufficient to resolve the structure of the flow 
and the internal structure of the shocked shell. The simulations 
were run for 0.1 Myr, after which the star has lost 10 Mq for 
M - 10 ^"' Moyr^' or 2 Mq for M = 2 x 10 ^^ Moyr^'. Eor 
Voo > 20 km s ' the simulations had reached a steady state at 
this time. The shocked shell has not reached a steady state for 
simulations with Voo = 15 kms^' (it is still accumulating mass), 
but this does not affect the ionized gas line emission. We checked 
this by running some simulations for 0.2 Myr and comparing the 
resulting emission properties. 


3.4. Structure of the photoionized wind 


The structure of the photoionized wind as a function of radius is 
shown in Eig.[3lfor simulation M5V15E. The ionization front is 


’ig-igf _ 

D-type ( denseTjKahn||1954 1 because the wind velocity satisfies 
Voo < 2ai, where Ui is the isothermal sound speed in photoion¬ 
ized gas. It has the classical structures of (from right to left) the 
ionization front, a dense shocked shell, and a shock in the neu¬ 
tral gas. Mackey et al. ( 2014| l assumed isothermal neutral gas, 
and so the shock was very compressive, whereas here the shock 
is adiabatic with a factor of 4 increase in density. Then there 
is a post-shock cooling region where the temperature decreases 
(from r ^ 0.043 - 0.0455 pc), and density increases to compen¬ 
sate (velocity also decreases to conserve mass). Einally there is 
a dense cooled shell bounded by the cooling layer and the ion¬ 
ization front. 

Simulations with Voo = 20 and 25 km s ' are shown in 
Figs. 0 and 1^ respectively. Eor Voo = 25 km s ' the wind is 
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Fig. 3. Radial structure of the wind in simulation M5V15E: we show 
the H number density fraction, gas radial velocity Vr, and tem¬ 

perature r on a logarithmic scale as a function of distance from the red 
supergiant. 



Position (pc) 

Fig. 4. As Fig.|^but for simulation M5V20E. 


too fast to sustain a shocked region and the ionization front is 
R-type ( rarefied ), characterised by weak density and velocity 
changes across the front ( Kahn|1954 i. The temperature increases 
smoothly in the wind from a few hundred to about 10 000 K. The 
ionization front is also much broader because there is no dense 
shell to absorb all of the photons. Simulation M5V20E (Fig. 
is a transition case on the boundary between D-type and R-type 
ionization fronts. There is a shock in the wind, but the shock can¬ 
not propagate far upstream from the ionization front and so there 
is too little time for the post-shock gas to cool. The separation be¬ 
tween the shock and the ionization front is calculated in |Macke^ 
et al.|(2014l and goes to zero as Voo ^ 2ai. In this case there is a 


shocked shell, but it cannot accumulate significant mass. 

Figs. [3||^ also plot which determines (along with tem¬ 
perature) the emissivity of the nebular emission lines. The peak 
value of Mg is always right at the ionization front, and drops by 
an order of magnitude between r - 0.044 pc and r = 0.065 pc. 
We therefore expect most of the line emission to arise from a 
relatively thin shell of ionized gas around the ionization front. 


3.5. Properties of the Ha and [NII] line emission 

We project the spherically symmetric wind onto the plane of the 
sky and trace rays through the wind for various impact param¬ 
eters, b, the radius of closest approach for each ray to the red 
supergiant. The emission and absorption coefficient for the Ha 
and [Nil] lines are calculated for every grid zone and tabulated. 
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Fig. 5. As Fig.j^but for simulation M5V25E. 


Then each ray is split up into line segments corresponding to the 
grid zones, and the radiative transfer equation is solved analyt¬ 
ically along each segment assuming a constant source function 
(see e.g. Rybicki & Fightman|1979| l, first inwards Xo r - b and 
then back outwards. To calculate velocity information we obtain 
the line-of-sight radial velocity of the gas from the gas dynamics, 
and the emission is then thermally broadened about this value us¬ 
ing the gas temperature. We use 501 velocity bins between -50.1 
and 50.1 km s ' to ensure that the lines are always well resolved 
in velocity. 

The radial velocity of peak emission is then extracted for 
each pixel, and the line brightness is obtained by summing the 
emission at all velocities for each pixel. The line ratio is the ra¬ 
tio of these line intensities. We convert from impact parameter, 
b, to an angular distance from W26, 9, by assuming a distance 
to Westerlund 1 of 3.55 kpc ( |Brandner et al.|[2008[ ), consistent 
with the distance derived more recently from eclipsing binary 
stars (3.7 ± 0.6 kpc) by Koumpia & Bonanos ( 2012| l. We also 
spatially smooth the intensity maps to a resolution of 1 arcsec at 
the distance of Westerlund 1, for easier comparison with W26. 

The results are plotted in Fig.j^for simulation M5V15F, with 
the observational data overplotted, as a function of distance from 
W26 measured in arcminutes. Fine intensities for [Nil] and Ha 
are plotted in panel (a), the line ratio in panel (b), and radial 
velocity of the line peaks in panel (c). The relative normalisation 
of the simulation to the data is chosen by setting the mean [NII] 
line intensity to unity between the two radii shown by vertical 
dotted lines in panel (a). Note again that the simulation results 
assume nitrogen is enhanced in the stellar wind by a factor of 
2.8 (see section]^. The spectral data close to W26, within about 
0.015 arcmin, are less reliable than at larger separations because 
of stellar continuum subtraction (see section]^. 

This simulation has a D-type ionization front with a dense 
static shell, so the ionization front is much denser than in the 
higher velocity simulations. This leads to a thinner region of 
peak emission because the newly ionized gas accelerates, and so 
its density (and hence emissivity) decreases more rapidly with 
radius than for a constant velocity wind. This simulation there¬ 
fore has the most limb-brightened emission, compared with sim¬ 
ulations with faster winds. 

The peak emission in panel (a) is at the same position as the 
observations by construction, because we constrained the ioniz¬ 
ing flux to locate the ionization front at this radius. The observed 
emission is much less peaked than the synthetic data for both 
lines, although this can partly be explained by the finite spatial 
resolution of the observations. The degree of limb brightening 
constrains the thickness of the shell of brightly emitting gas, with 
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Position (arcmin) 

Fig. 6. The Ha and [Nil] spectral lines from the simulation M5V15E, 
with M = 2 X 10^^ Moyr-', Voo = 15 hms^', and = 1.70 x 
10*° cm“^ s“*, after 0.1 Myr of evolution. The dashed lines show the 
simulation, solid lines smoothed to 1 arcsec resolution at the distance 
of Westerlund 1, and the dotted lines with points show the observations 
to the ESE of W26. The panels show (a) line brightness, (b) line ratio 
[NII]/Tlti', and (c) radial velocity of the peak emission, all as a func¬ 
tion of distance from the star. Vertical dotted lines in panel (a) show 
the range of radii used to set the normalisation of the simulated and 
observed [NII] line brightness. 


a thinner shell being more limb-brightened. The data thus appear 
to imply a thicker shell of brightly-emitting gas compared with 
what we find from this simulation. 

The [Nil] brightness decreases more rapidly with distance 
from W26 than the synthetic data, but this is less clear with the 
Ha brightness because the ratio between maximum and mini¬ 
mum observed brightness is not as large as for [Nil]. The ob¬ 
served diffuse Ha emission at large radii is a constant component 
from Westerlund 1 and not related to the W26 nebula (see Fig.|^, 
even though its brightness is similar to our predictions at sepa¬ 
ration 0.07 -0.1 arcmin. Subtraction of this constant component 
complicates any estimate of the rate of decrease at these larger 
radii, but it may be that the observed Ha also decreases more 
rapidly with radius than we predict. The rate that the brightness 
decreases with distance from the peak constrains the density of 
the flow as a function of radius if the gas is isothermal (assuming 
spherical expansion). The data then indicate that either nitrogen 
begins to be doubly ionized further from the star, or the emitting 
gas density decreases much more rapidly than our models pre¬ 




Fig. 7. The Ho- and [Nil] spectral lines from the simulation M5V25E, 
again after 0.1 Myr of evolution. The lines and symbols are the same as 
in Fig.|^ 


diet. This could arise if the wind of W26 is hydrodynamically 
confined by much hotter intracluster gas (see Sect.|^. 


The brightness ratio, [NII]/Ha, is plotted in panel (b) and 
shows that within the bright part of the nebula the synthetic and 
observed data have a broadly constant line ratio that is 2 for the 
observations and a; 2.5 for the simulations. If the wind were not 
enriched in nitrogen, the simulations would predict a line ratio 
of ai 0.9 in this region, inconsistent with the results. The obser¬ 
vations therefore indicate that nitrogen is enhanced by a factor 
of ~ 2 - 3 in the nebula. If we consider the nitrogen abundance 
to be a free parameter, we obtain good agreement between sim¬ 
ulations and observations for f{N) - 1.60 x 10“"^, or 2.35 times 
the solar abundance ( Asplund et al.|2009 1 . This enhancement is 
strong evidence that the nebula contains gas processed by the 
CNO cycle, and hence that it consists of mass lost by W26, at 
least out to a distance of0.06 arcmin. 


The line ratio decreases in the simulations from the ioniza¬ 
tion front (0 0.045 arcmin) outwards, whereas the observed 

line ratio stays constant to the edge of the bright nebula and then 
decreases sharply from « 1.5 to < 0.5 between 0 = 0.06 and 
0 - 0.08 arcmin. The synthetic data never decrease to such a 
small line ratio, so again we suggest for this low-emission gas 
that either nitrogen is more highly ionized or the gas is not en¬ 
riched in nitrogen (i.e. it is not part of the wind of W26). What¬ 
ever the explanation, clearly there is something happening in the 
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Fig. 8. Comparison between wind velocity, v„, and the blueshift of 
the Ha and [Nil] spectral lines for lines of sight close to the star. 
The boundary between D-type and R-type ionization fronts is at v„ « 
22 km s“*. The red lines are for Ha and the blue for [NII]. 


region between 6 - 0.06 and 0 — 0.08 arcmin that is not seen in 
the synthetic data. 

This is further shown by looking at the radial velocity of the 
line peak as a function of position, shown in panel (c). The in¬ 
ner part of the nebula has radial velocity Viad ~ -15 kms“\ 
similar to the simulated lines. The observations then show that 
the brightest emission at 0.035-0.055 arcmin is also the most 
blueshifted, whereas the simulations show that the radial veloc¬ 
ity of the brightest emission is zero. This is the biggest dis- 
crepency between our model and the data: both the line ratio 
and the radial velocity stay approximately constant throughout 
the bright nebula, while the simulations show a change at the 
brightest part of the nebula. The decrease in observed line ratio 
is spatially correlated with the decrease in radial velocity. 

Fig. [7] shows the same plots for simulation M5V25E, with 
a faster wind such that the ionization front is R-type, with only 
weak density and velocity changes in the wind (see Fig.|^. The 
agreement of the line brightness with observations is better than 
for M5V15E, in that the line emission is not so strongly peaked 
at the ionization front, but the decrease in intensity at 0 > 0.05 
arcmin still has the wrong shape for [Nil]. The line ratio also 
agrees well for 6 < 0.06 arcmin. We see the same discrepen- 
cies for this simulation as for M5V15E, however. The simulation 
predicts that the line goes to zero radial velocity at the emission 
peak, whereas the data show the maximum radial velocity at the 
peak. The same holds true for the line ratio, which in the ob¬ 
servational data remains large as long as the radial velocity is 
signihcantly different from zero. 

The situation is somehat different on the opposite side of 
W26, Westnorthwest (WNW) in the direction of the extreme 
sgB[e] star W9. Here the nebula is conhned closer to W26, pos¬ 
sibly because of the proximity of W9. The observed data shows 
no limb-brightened peak in the nebular emission, rather just a 
shoulder of emission between 0.02 and 0.03 arcmin. The results 
of a comparison between simulation and observation are very 
similar to what we have already discussed for the ESE side of 
the nebula, and are shown in Appendix [B] 


3.5.1 . Wind velocity as traced by Ha and [N II] lines 

Eig. plots the blueshift of the two lines as a function of wind 
velocity for the four sets of simulations (two values of M and two 
sets of ionizing fluxes). Eor slow winds with D-type ionization 
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fronts (voo ^ 22 kms^') the ionized gas is accelerating outwards 
from the outer edge of a static shell, and so the line velocity is 
independent of wind velocity. In this velocity range the Ha line 
is blueshifted by about 2 kms^' more than the [Nil] line for 
all simulations. This is a temperature effect: Eq.[^shows that the 
two lines will be produced from different depths in the ionization 
front (see e.g. [Hester et al.|1996[ [Henney et al.|2005| . We also 
hnd that the simulations with denser winds (M - 10“'^ Mq yr ') 
are blueshifted by about 2 km s ' compared with simulations 
that have less dense winds. This is a density effect because within 
the ionization front, the heating rate scales linearly with density, 
whereas the cooling rate is quadratic for collisional coolants, 
implying that the temperature structure within the front should 
be density-dependent. The ionization front thickness is inversely 
proportional to density, and so its internal structure is also ex¬ 
pected to depend on M. We hnd that for a given temperature in 
the photoevaporation how, the denser winds have an expansion 
velocity about 2 km s ' larger than the less dense winds. 

This result for D-type ionization fronts likely depends 
weakly on the temperature of the ionizing radiation held: harder 
radiation helds create hotter ionization fronts, faster photoe¬ 
vaporation hows, and hence the lines will probably be more 
blueshifted. The metallicity of the wind is also important, be¬ 
cause lower metallicity gas has a higher equilibrium tempera¬ 
ture, which will also drive faster photoevaporation hows. 

Eor faster winds (voo ^ 22 km s ') the ionization front be¬ 
comes R-type, and this has very little effect on the wind dy¬ 
namics except for a slight deceleration. The ionized gas exiting 
the ionization front has a velocity given by (jMihalas & Mihalas 
[T9841 ch. 106) 


1 H- 


2\ 


(3) 


where the term {ajvoof' < 0.25 for R-type ionization fronts. This 
explains why the line velocity is always less blueshifted than Voo, 
and also why both lines approach Voo for larger Voo. Again, [NII] 
is less blueshifted than Ha because it forms deeper in the ion¬ 
ization front before thermal re-acceleration of the wind can take 
place. We conclude that these lines are useful but imperfect trac¬ 
ers of the wind velocity for red supergiants, with their accuracy 
increasing as wind velocity increases. 


3.5.2. Limb brightening of the Ha and [N II] lines 

Figs. 0 and|7 show that the slow wind simulations have more 
limb-brightened emission than faster winds. In Eig. [^a) this is 
quantihed with plots of the ratio, tlb, of the peak to the central 
intensity for all 16 simulations and for both spectral lines. This 
conhrms the trend that the two slow-wind simulations (Voo = 15 
and 20 km s ') with D-type ionization fronts have stronger limb 
brightening than the two fast-wind simulations (voo = 25 and 
30 km s ') with R-type ionization fronts. The slow-wind simu¬ 
lations have tlb ~ 3.5 - 5.5 for Ha and tlb ~ 4.5 - 6.5 for 
[Nil]. The fast-wind simulations have tlb ~ 3.3 - 4.3 for Ha 
and tlb ~ 4 - 5.1 for [Nil]. This arises, as noted above, because 
the emitting region is thinner for the slow-wind simulations, and 
limb brightening is stronger the thinner the emitting shell. [NII] 
emission arises from a thinner shell of emitting gas within the 
ionization front than Ha emission, and so it has larger tlb than 
Ha for the same reason. 

The simulations with denser winds have smaller tlb, which 
arises because of self-absorption. Lines of sight towards the star 
have very little absorption (from the front side; the back side is 
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Fig. 9. Comparison between wind velocity, Voo, and the degree of limb 
brightening of the Ha and [NII] spectral lines, defined as the ratio of the 
peak brightness to that at lines of sight near the star (r = 0). The bound¬ 
ary between D-type and R-type ionization fronts is at v^c ~ 22 kms“*. 
The red lines are for Her and the blue for [Nil]. Panel (a) shows the raw 
simulation data, and panel (b) when smoothed to a spatial resolution of 
1 arcsec. The lines in panel (b) follow the legend in panel (a). 


invisible) because the wind density decreases outwards. Lines 
of sight at the limb pass through much more material, however, 
and the opacity is significant for the simulations with the densest 
winds. 

Fig. I^b) shows, however, that the limb brightening is much 
decreased when the data are smoothed to a resolution of 1 arc- 
sec, with most values of tlb lying between 2 and 3.5. We still 
see the trend that denser winds have less limb-brightening, but 
the difference between Ha and [Nil] is almost completely gone. 
Also, the trend of decreasing tlb with increasing Voo is no longer 
present, and is even reversed for simulations with dense winds. 
The smoothed values are still somewhat larger than the observed 
ratio for W26, z-lb ~ L5. 


4. Discussion 

4. 1. Constraining wind properties using photoionized winds 

Photoionized red supergiant winds can be extremely bright, be¬ 
ing very strongly irradiated by nearby stars. The peak bright¬ 
ness of the simulated nebulae in Ha (smoothed to a spatial res¬ 
olution of 1 arcsec for the distance of Westerlund 1) ranges 
from, /max ~ 10“'^ ergem^^ s“' arcsec^^ for M5V30E to /max ~ 
10“^' ergem^^ s“* arcsec^^ for M4V20W. This is, for example, 
200 - 20000 times brighter than the photoionized bow shock 
around the red supergiant IRC -10414 ( [Meyer et al.|2014al l, and 


=» 1-100 times brighter than the Ha emission from the bright¬ 
est of the Eagle Nebula pillars ( Hester et al.|1996l l. The overall 
brightness scales approximately with the square of the electron 
density (hence the wind density). This can potentially provide a 
good constraint on the mass-loss rate of red supergiants, in the 
same way that radio data gives the emission measure and hence 
electron density (with assumptions about geometry). 

Using the brightness estimate from jWright et al.] ( |2014| l for 
the circumstellar nebula around W26, and assuming the same 
extinction as they find for the nearby triangular nebula (8.4 mag 
at Ha), the mean intrinsic brightness of the circumstellar nebula 
is /(Ha) ~ 6 X 10 ergem^^ s ' arcsec^^. This puts the nebula 
at the upper end of the emission range of our models, and may 
be evidence that the mass-loss rate is closer to our upper value 
than our lower value. Our simulations have similar electron den¬ 
sities (few xlO^^ cm“^, see Eigs. 3]|5 i to the radio estimates from 
[Dougherty et al.| ( |2010| l, so this agreement is not surprising. This 
holds even though our model does not match the data very well, 
because the mass-loss rate is the primary parameter determining 
the gas density at a given distance from the star. 

Winds that are slow enough to allow D-type ionization fronts 
(< 22 kms * in these models) have blueshifted lines for which 
the velocity is determined by the ionization front and not by the 
wind. We find line blueshifts of 16 - 19 km s ', although this de¬ 
pends somewhat on wind density (and probably also on the radi¬ 
ation temperature). Eor faster winds an R-type ionization forms, 
and this has only weak dynamical effects on the wind. In this 
regime, nebular lines become good tracers of the wind veloc¬ 
ity, although the line velocity is always an underestimate. This 
is potentially very useful for estimating wind terminal velocities 
because the best estimates (usually from masers) come from the 
wind acceleration region itself and require models for the turbu¬ 
lent velocities in the wind (e.g. for NML Cyg, [Zubko et al.[2004[ 
Nagayama et al. 2008[l. Eok et al. (2012[) detected maser emis¬ 


sion from the wind of W26, but were not able to measure the 
wind velocity with their data. Advances in instrumental capabil¬ 
ities will soon enable measurements of the wind velocity using 
similar methods as for NML Cyg, so an independent constraint 
from nebular emission at larger radii is a useful check on these 
models. 


4.2. Comparison with \N26 

We have shown that the there is some agreement between ob¬ 
servations and the predictions of the spherically symmetric pho¬ 
toionized wind described in [Mackey et al.[ ( [2014| l, but that there 
are also important differences that seem difficult to reconcile 
with any spherically symmetric wind model. The whole pho¬ 
toionized nebula seems to be significantly blueshifted wherever 
there is significant emission, favouring an interpretation where 
W26 is on the near side of Westerlund 1 and its wind is being 
swept towards us by some process (either radiative or hydrody- 
namical) associated with winds, radiation, or supemovae from 
the star cluster. 

This seems more similar to the one-sided photoionization 
of the wind of NML Cyg (by the Cyg OB2 association; [Mor-[ 
ris & Jura|1983] l than to the diffuse ionizing radiation field that 
is invoked to explain the photoionization of Betelgeuse’s wind 
( [Mackey et al.|2014| . If the [Mackey et al.[ ( [2014[ ) model is cor¬ 
rect, then Betelgeuse is irradiated by an isotropic ionizing radi¬ 
ation field, creating a spherically symmetric ionization front in 
the star’s wind and producing a static shell. Betelgeuse is in this 
relatively quiet environment because it is a runaway star that is 
now far from its place of birth. NML Cyg is very different, being 
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near a massive and young OB association, and so experiencing a 
very anisotropic radiation field. W26 experiences an even more 
extreme environment because it is still within (at least in projec¬ 
tion) the most massive star cluster in our Galaxy, that contains 
many WR and O-type stars. It follows then, that we should not 
be surprised if W26 were to experience an asymmetric radiation 
field, and if a multi-dimensional extension of the [Mackey et al.| 
( 2014| l model were required to explain the data. We discuss this 
and other alternative explanations in the following subsections. 


shocked wind of the WR and O stars in Westerlund 1, by the 
numerous supernovae that have already occurred, or a combina¬ 
tion of both. These authors estimate the H number density and 
temperature of the X-ray emitting plasma to be hh = 0.5 cm~^ 
and kT - 3.7 keV (T - 4.3 x 10^ K), respectively. For a fully 
ionized plasma containing 1/11 helium by number, this gives a 
thermal pressure of po - 2.2nukT - 6.5 x 10“® dynecm^^. Bal¬ 
ancing this against the ram pressure of the wind from W26 gives 
the radius of the wind termination shock, Rts, of 


4.3. Bow shock 


If the nebula is a bow shock then it must be oriented away from 
us, i.e., W26 is embedded in an ISM with a bulk flow directed 
approximately towards us. The nebula is blueshifted by about 
25 km s ' relative to W26, which represents the velocity of the 
wind material that has been entrained in the bulk flow by the 
bow shock. This radial velocity then measures the relative ve¬ 
locity between the star and the ISM. The standoff distance, Rgo, 
is then less than or equal to the smallest radial extent of the neb¬ 
ula, which we estimate to be 0.037 pc, and is given by 


^so - 



(4) 


where po is the ISM density and vq = 25 km s“' the flow ve¬ 
locity of the ISM past W26. The wind is photoionized, so we 
assume this accele rates the ionized part of the wind to at least 
Voo = 35 kms^' (Meyer et al. 2014aI. We then use our up¬ 
per and lower estimates of M for W26 to estimate the value 
of Po that is required to produce the observed nebula. We find 
Po = 2.15x10“^" gcm“^ (or no = 10"'cm"^, where no is the num¬ 
ber density of H and He nuclei) for M - 10 Mq yr ', or 5 times 
smaller for M = 2 x 10“^ Mq yr“'. This seems an unreasonably 
large ISM density for the intracluster medium of Westerlund 1 
(X-ray emission suggests no ~ 0.5 cm“^; see below), arguing 
against a straightforward bow shock interpretation for the neb¬ 
ula. For reference to the following discussion, the ram/thermal 
pressure required to confine the wind to Rso is p = 1.3 x 10"^ 
dyne cm“^ if M - 10“"' Mq yr“', or p = 2.7 x 10“^ dyne cm“^ if 
M = 2 X 10“^ Mq yr“', both calculated for Voo = 35 km s“'. 

To further test this model, we relax the assumption that the 
shocked wind of W26 is accelerated to the full velocity of the 
bulk flow, which might be possible if there is little mixing at 
the contact discontinuity. This then allows a much faster bulk 
flow with a lower density, but to decrease the flow density to 
no ~ 1 cm“^ we require a gas velocity, vo > 10^ kms“'. This 
sort of gas flow in the intracluster medium is only likely in the 
immediate vicinity of an extreme driver of momentum, e.g. a 
massive star wind or a supernova. In comparison to other mas¬ 
sive star clusters, Povich et al.j ( |2008[ ) estimate ram pressures 
of poVq ~ (0.2 - 5) X 10 ® dyne cm“"' for bow shocks around O 
stars in the M17 and RCW 49 massive-star-forming regions. This 
shows that the nebula around W26 is probably even more ex¬ 
treme than any of the bow shocks studied by Povich et al. ( 2008| l, 
if the bow shock interpretation is correct, suggesting again that 
perhaps we need an extreme source of momentum very nearby, 
rather than a cluster-scale gas outflow. 


4.4. Thermal pressure confinement 

Thermal pressure confinement is a promising mechanism, be¬ 
cause X-ray observations have revealed a hot, diffuse, intraclus¬ 
ter medium ( [Kavanagh et al.|[2011| l. This is produced by the 
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Ris - 0.17 pc 


M 

yj 10-"'Moyr-' 


Vo< 


35 kms 


(5) 


For our two mass-loss rates we obtain Rts = 0.17 pc (M = 
10“"' Moyr“') and Rts = 0.075 pc (M = 2 x 10“^ Mgyr''), 
both of which are larger than the observed extent of the neb¬ 
ula (0.037 pc). For the lower mass-loss rate, however, Rts is not 
much larger than the observed nebula, and a factor of 4 increase 
in the external pressure would bring the numbers into agreement 
(or a 4x decrease in M to 5 x 10“® Moyr“', although this is 
unlikely for such a massive red supergiant). Some asymmetric 
pressure is still required, for example a subsonic bulk flow (e.g. 
[Mackey et al.|2015] l, to create the asymmetric nebula, so purely 
thermal-pressure confinement is not a viable explanation. 


4.5. Wind-wind interaction 


The star W9 is 0.23 arcmin from W26 in projection, or 0.24 pc 
at the distance of Westerlund 1. It is the strongest radio emitter 
in Westerlund 1 and is surrounded by a dense wind with M 
3.3 X 10“"'Mq yr“' and Voo ~200kms“' (Dougherty et al. 2010 


model 3). If the winds of the two stars collide, then the wind of 
W26 will be shocked and decelerated at a distance, R^^i from 
W26 that is given by 


Rwc ~ 


1 - Vc 

1 - c 


(6) 


where Rq is the separation of the two stars and c is the ratio of 
the wind momenta: 


(MVoc) 

V /W26 

(mvJi 

V /W9 


(7) 


For the W9 wind properties above, and using M - 10“"' Mq yr“' 
and Voo = 35 kms ' for W26, we find = 0.19Ro- For W26 
with M - 2x 10“^ Mq yr“', this reduces to Rwc = 0.093Ro- The 
observed radius of the nebula is 0.037 pc, implying Rq = 0.20 
pc for the larger M estimate or Rq - 0.40 pc for the smaller M. 
This is again comparable to the observed separation between the 
two stars, and so should be considered as a plausible alternative 
scenario to a simple photoionized wind. 

Considering the pressure of the X-ray plasma quoted above 
(6.5 X 10“^ dyne cm“^), the wind of W9 would have a termina¬ 
tion shock at Rts = 0.73 pc. The freely expanding wind region 
would have a diameter of 1.4 arcmin at the distance of Wester¬ 
lund 1, covering a large fraction of the core of Westerlund 1. This 
means that W26 and many other stars in the cluster may well be 
embedded in the freely-expanding wind of W9, and could be ex¬ 
periencing a wind-wind collision. 
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4.6. Photoionized wind 

We have shown that a spherically symmetric model for a pho¬ 
toionized wind has severe problems explaining the observations, 
but perhaps an asymmetric model could work. This was pro¬ 
posed by [Morris & Jura ( |1983| l for NML Cyg, to explain the 
arc-shaped H ii region on one side of the star. In this model there 
is one dominant radiation source so that the wind is photoionized 
mainly from one direction and the ionization front has a shape 
similar to a bow shock (for mathematical details, see [Morris &| 
[Jura[[198^ . If the ionization front is D-type then the wind can 
be decelerated into a one-sided shell, and the asymmetry of the 
shell leads to non-radial flows away from the radiation source 
because of the rocket effect ( Oort & Spitzer[1955] l. This requires 
multi-dimensional simulations to model quantitatively, beyond 
the scope of this paper, but we intend to investigate this in more 
detail in future work. 

The EUV fluxes in Table show what is required to ionize 
the wind of W26, and simple estimates show that nearby stars 
can provide these fluxes. The Wolf-Rayet star WR77o (WR B) 

( Crowther et al. 12006 1 is about 0.5 pc from W26 (projected dis¬ 
tance). According to Crowtherj ( [2007[ l, it should have an ioniz¬ 


ing photon luminosity, Qo - Iw^^s '. For a distance of 0.5 
pc, the EUV flux at W26 is Fy = 8.4 x lO" cm^^s^', com¬ 
parable to the largest fluxes used in the spherically symmet¬ 
ric simulations. Its wind, by contrast, is probably not strong 
enough to shape the nebula, because Crowtherj ( [2007[ ) estimate 
M - 1.6 X 10“^ Moyr^' and Voo = 1300 kms^\ too weak to 
confine W26’s wind according to Equation (|^. 

Much closer to W26, the 09 lab star W25 h as an ionizing 
photon luminosity of Qo - (lO"**'^ - 10"^® °) s"^ (Martins et al. 


|2005| l, and is only 0.09 pc from W26 in projection. This corre 
sponds to an ionizing flux of Fy - (6.4 - 10.0) x lO'^ cm“^ s“\ 
and so could ionize the wind of W26 even deeper than the ob¬ 
served bright nebula, if the physical separation of the two stars 
is comparable to the projected separation. Another candidate is 
W9, discussed above in the context of wind-wind collision, but 
we do not have any reliable estimates of its ionizing photon lu¬ 
minosity. 

4.7. Hot companion star? 

Westerlund 1 was observed with Chandra in X-rays and analysis 
of point sources was presented by Skinner et al. ( 2006[ l. W9 was 
clearly detected, which they argue is because W9 is a colliding- 
wind binary, but W26 (referred to as Ara A) was not detected 
although it was observed. If W26 was a binary system with a 
neutron star (or black hole) companion that provided the ioniz¬ 
ing photons to produce the nebula, then it would also be bright 
in X-rays, and this is not observed. W26 is also too young to 
host a white dwarf companion, so we can exclude that a com¬ 
pact companion is ionizing the wind. It is also not a photometric 
variable star ( [Bonanos[2007[ ) and, although it has spectral varia¬ 
tions, these are consistent with other massive red supergiants and 
not indicative of duplicity (Clark et al.|2010|l. 


4.8. Outiook 

We see that there are many radiation sources that should pho- 
toionize the wind of W26 to at least the radius of the observed 
nebula, and potentially much closer to the star, although it is not 
yet clear if one-sided photoionization could produce the veloc¬ 
ity signatures of the observed nebula. We have also shown that 
the measured X-ray plasma pressure could contribute to confin¬ 


ing the wind. The large temperature difference between the pho¬ 
toionized red supergiant wind and the X-ray plasma of the in¬ 
tracluster medium may also set up a strong thermal conduction 
front at the contact discontinuity between wind and ISM, which 
has been shown to strongly modify the properties of O star bow 
shocks ( Comeron & Kaper|1998[ Meyer et al. [20145) 1. We sug¬ 
gest, however, that a wind-wind collision between the winds of 
W9 and W26 is the most obvious scenario that could produce 
the observed nebula. Naively, we expect that a bow shock model 
(possibly from a wind-wind collision with the wind of W9) is 
the only one that can explain why of the bright emission is 
blueshifted. 

It would be very useful to have higher resolution spectra, 
to reduce the uncertainties in the radial velocity of the spec¬ 
tral lines, and to search for substracture in the lines (e.g. red- 
shifted absorption). Better spatial coverage of the nebula is also 
required. In particular spectra of the triangular nebula discov¬ 
ered by Wright et al. ( 2014[ l and the emission connecting it with 
the W26 nebula could reveal velocity gradients with distance 
from W26. This would show whether the gas is being acceler¬ 
ated away from W26, as expected for gas that is entrained by 
a wind. This would also provide more conclusive evidence for 
whether W26 is on the near or far side of Westerlund 1: if on the 
near side then the gas should be more blueshifted the further it is 
from W26. 


5. Conclusions 

We have presented the first predictions for nebular emission from 
photoionized winds of red supergiants, using spherically sym¬ 
metric, radiation-hydrodymamical simulations and postprocess¬ 
ing radiative transfer. Although the results have been compared 
only to W26, they also have application to photoionized winds 
around other red supergiants. The simulations improve on the 
previous work of ( [Mackey et al.[2OT4) l in that the gas temper¬ 
ature is calculated self-consistently from mechanical and radia¬ 
tive heating/cooling, allowing us to make realistic predictions of 
temperature-sensitive nebular emission from ionized gas. 

We produce position-velocity datasets for Ha and [Nil] 
spectral lines from the simulations, and study the properties 
of this emission. The peak brightness of the simulated nebu¬ 
lae is 10“'^ - 10“" ergcm^^ s“'arcsec^^, for mass-loss rates. 


M = (0.2-l)x 10“^ Moyr“‘, and wind velocities, Voo = (15-30) 
km s '. Simulations of slow winds that are decelerated into a 
dense shell show strongly limb-brightened line emission, and ra¬ 
dial velocity of spectral lines is independent of the wind speed. 
Faster winds (> 22 km s ') that cannot form a dense shell have 
less limb-brightening and the line radial velocity is a good tracer 
of the wind speed. All of the simulations have the most blue¬ 
shifted emission at the position of the star, and the radial velocity 
of the lines goes to zero at the region of brightest emission (the 
limb), as is typical for emission from an expanding shell. 

If nitrogen is not enriched in the wind, then the predicted 
line ratio of [Nil] to Ha is about 0.9. Nitrogen is enhanced in 
red supergiant winds of a factor of about 3 (cf. [Brott et al.[20lT| l, 
and for this we obtain line ratios of about 2.5 - 2.8. 

We have compared spectra of the Ha and [Nil] emission 
lines from the nebula around the red supergiant W26 in Wester¬ 
lund 1 ( Wright et al.[[2014) l with our simulations to investigate 
the origin and physical properties of the nebula. We also obtain 
a temperature for the star Teff = 3600 K from the same spectra. 
The brightness distribution of the emission and the line ratio of 
[NII] to Ha show similarities to the photoionization model, but 
there are discrepencies for the radial velocity of the emission. 
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All of the observed bright nebular emission is blueshifted with 
respect to W26 by 15-25 kms“\ and the brightest emission is 
the most blueshifted, neither of which can be reproduced by a 
spherically symmetric model of a photoionized wind. From the 
observed line ratio (> 2.25), we can conclude that the nebula 
is photoionized, however, and that it is enriched in nitrogen by 
a factor of about 2.35. This shows that the nebula is produced 
by the wind of W26, although the discrepencies with our model 
suggest that the outer edge of the nebula may be hydrodynami- 
cally confined rather than freely expanding. 

We propose that a wind-wind collision scenario, where the 
wind of W26 is swept towards the observer by the stronger wind 
of the nearby star W9, may match the data better. Our results 
suggest that the circumstellar medium may have more struc¬ 
ture even closer to the star than the observed nebula, because it 
does not appear as though we have detected the ionization front 
and its associated shell in the current observations. Multidimen¬ 
sional models and/or simulations, and spectra with higher resolu¬ 
tion and better spatial coverage are necessary to progress further 
in understanding the circumstellar nebula around W26. Further 
study of this nebula will contribute to our understanding of what 
happens to the mass lost by red supergiants in star clusters. It is 
important to know whether it is rapidly heated and expelled from 
the cluster, or whether it can remain within the cluster for an ex¬ 
tended period of time, remain cool and dense, and potentially be 
available for secondary star formation in the cluster. 
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Appendix A: Results for the ESE side of the nebula 


Here we plot the line intensity, radial velocity, and line ra¬ 
tio for the simulations that are not discussed in the text. 


Eigs. A.l and 

A.2 

show M5V20E and M5V30E, respec- 

tively, and Eigs. 

A.3 

A.6 show results for M4V15E, M4V20E, 


M4V25E, and M4V30E. 
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Position (arcmin) 

Fig. A.l. The Via and [Nil] spectral lines from the simulation 
M5V20E, with M = 2x 10^^ Moyr-‘, v„ = 20 hms^', and Fy = 
2.78 X 10*® cm“^ s“*, after 0.1 Myr of evolution. The solid lines show 
the simulation, and the dotted lines with points show the observations 
to the ESE of W26. The panels show (a) line brightness, (b) line ratio 
[N Ili/Hfr, and (c) radial velocity of the peak emission, all as a function 
of distance from the star. 


Appendix B: Results for the WNW side of the 
nebula 

Here we plot the line intensity, radial velocity, and line ratio for 
the simulations labelled “W” in Table [T] These model the neb¬ 
ula Westnorthwest (WNW) from W26, in the direction of the ex¬ 
treme sgB[e] star W9. As noted in the text, the nebula is conhned 
closer to W26 on this side, possibly because of the proximity of 
W9. The “W” simulations have larger photon fluxes than the cor¬ 
responding “E” simulations, so that the ionized region penetrates 
deeper into the wind of the red supergiant. The results of simu¬ 
lations are overplotted on the observed data in Fig.|B.l|for sim¬ 
ulation M5V15W, and similarly in Figs.|B.2HB.4|for M5V20W, 
M5V25W, and M5V30W. Figs.|B.5p.8|show results for the sim¬ 
ulations M4V15W, M4V20W, M4V25W, and M4V30W, which 
have a denser wind because of a larger mass-loss rate. 

The observed data shows no limb-brightened peak in the 
nebular emission, rather just a shoulder of emission between 
0.02 and 0.03 arcmin. If we discard the uncertain two points 
closest to zero (at W26) then the emission decreases monotoni- 
cally with distance from W26, in strong constrast to the predic¬ 
tions of simulations. This makes the discrepency with the limb- 
brightened emission from M5V15W (Fig. |B.l^) even larger for 




Position (arcmin) 


Fig. A.2. The Ha and [Nil] spectral lines from the simulation 
M5V30E, with M = 2x 10^^ Moyr^*, v,*, = 30 kms-‘, and Fy = 
1.26 X 10'® cm“^ s“*, after 0.1 Myr of evolution. The lines and symbols 
are the same as in Fig.|A.l[ 


this side of W26. The observed line ratio is slightly lower on this 
side of W26, but again is relatively constant across the nebula, 
decreasing only outside the bright part of the nebula {9 > 0.04 
arcmin). Additionally, the radial velocity of the lines (Fig. B) 
appears too low at all radii, but particularly for 6 > 0.03 arcmin. 

Simulation M5V25W (Fig. m has better agreement with 
observations, for the same reasons as M5V25E: the ioniza¬ 
tion front emission is less peaked, the radial velocity is more 
blueshifted in agreement with observations, and the line ratio 
has reasonable agreement with observations. Again, however, 
the crucial problem is that the observed nebular emission is all 
blueshifted by about 25 km s ', whereas in the simulation the 
brightest emission is neither blue- nor redshifted. 
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Position (arcmin) 

Fig. A.3. The Hor and [Nil] spectral lines from the simulation 
M4V15E, with M = lO^^^ Moyr-', Voo = 15 kms-‘, and Fy = 
3.67 X 10" cm^^s^*, after 0.1 Myr of evolution. The lines and sym¬ 
bols are the same as in Fig.|A.l| 



Position (arcmin) 

Fig. A.4. The Ha and [Nil] spectral lines from the simulation 
M4V20E, with M = lO^^^ Moyr-\ v<„ = 20 kms-‘, and Fy = 
5.05 X 10" cm^^s^*, after 0.1 Myr of evolution. The lines and sym¬ 
bols are the same as in Fig.|A. 1| 
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Position (arcmin) 

Fig. A.5. The Ho: and [Nil] spectral lines from the simulation 
M4V25E, with M = lO^^^ Moyr-', Voo = 25 kms-‘, and Fy = 
3.10 X 10" cm^^s^*, after 0.1 Myr of evolution. The lines and sym¬ 
bols are the same as in Fig.|A.l| 



Fig. A.6. The Ha and [Nil] spectral lines from the simulation 
M4V30E, with M = lO^^^ Moyr-\ v<„ = 30 kms-‘, and Fy = 
2.15 X 10" cm^^s^*, after 0.1 Myr of evolution. The lines and sym¬ 
bols are the same as in Fig.|A. 1| 
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Fig. B.l. The Hot and [Nil] spectral lines from the simulation 
M5V15W, after 0.1 Myr of evolution. The solid lines show the simula¬ 
tion, and the dotted lines with points show the observations to the WNW 
of W26. The panels show (a) line brightness, (b) line ratio [NII]/Ha, 
and (c) radial velocity of the peak emission, all as a function of distance 
from the star. Vertical dotted lines in panel (a) show the range of radii 
used to set the normalisation of the simulated and observed [NII] line 
brightness. 
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Fig. B.2. The Hot and [Nil] spectral lines from the simulation 
M5V20W, with M = 2x 10^^ Moyr-‘, v™ = 20 kms-‘, and Fy = 
1.29 X 10'* cm“^ s“*, after 0.1 Myr of evolution. The lines and symbols 
are the same as in Fig.|B.l[ 
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Fig. B.3. The Ha and [Nil] spectral lines from the simulation 
M5V25W, with M = 2 x 10^^ Moyr-', Vcx. = 25 hms^', and = 
8.29 X 10*** cm^^s^', again after 0.1 Myr of evolution. The lines and 
symbols are the same as in Fig.|B.l| 
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Fig. B.4. The Ha and [Nil] spectral lines from the simulation 
M5V30W, with M = 2x 10^^ Moyr-\ v™ = 30 kms-‘, and Fy = 
5.83 X 10'° cm“^ s^*, after 0.1 Myr of evolution. The lines and symbols 
are the same as in Fig.|B.l[ 
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Fig. B.5. The Hot and [Nil] spectral lines from the simulation 
M4V15W, with M = lO^"* Moyr-‘, v,,,, = 15 kms-‘, and = 
1.70 X 10'^ cm“^ s“*, after 0.1 Myr of evolution. The lines and sym¬ 
bols are the same as in Fig.|B.l| 
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Fig. B.6. The Hot and [Nil] spectral lines from the simulation 
M4V20W, with M = lO^^^ Moyr-‘, Voo = 20 kms-‘, and Fy = 
2.34 X 10*^ cm^^s^*, after 0.1 Myr of evolution. The lines and sym¬ 
bols are the same as in Fig.|B.l| 
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Fig. B.7. The Ha and [Nil] spectral lines from the simulation 
M4V25W, with M = lO^"* Moyr-‘, v,,,, = 25 kms-‘, and = 
1.44 X 10'^ cm“^ s“*, after 0.1 Myr of evolution. The lines and sym¬ 
bols are the same as in Fig.|B.l| 



Fig. B.8. The Ha and [Nil] spectral lines from the simulation 
M4V30W, with M = lO^^^ Moyr-‘, Voo = 30 kms-‘, and Fy = 
9.95 X 10** cm^^s^*, after 0.1 Myr of evolution. The lines and sym¬ 
bols are the same as in Fig.|B.l| 
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